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Abstract. In charge-stabilized colloidal suspensions, highly charged macroions, dressed 
by strongly correlated counterions, carry an effective charge that can be substantially 
reduced (renormalized) from the bare charge. Interactions between dressed macroions 
are screened by weakly correlated counterions and salt ions. Thermodynamic and 
structural properties of colloidal suspensions depend sensitively on the magnitudes 
of both the effective charge and the effective screening constant. Combining a charge 
renormalization theory of effective electrostatic interactions with Monte Carlo simula- 
tions of a one-component model, we compute osmotic pressures and pair distribution 
fxmctions of deionized colloidal suspensions. This computationally practical approach 
yields close agreement with corresponding results from large-scale simulations of the 
primitive model up to modest electrostatic coupling strengths. 

PACS (2006): 82.70.Dd, 83.70.Hq, 05.20.Jj, 05.70.-a 



1 Introduction 

Charge-stabilized colloidal suspensions [1-3] exhibit rich thermodynamic phase behav- 
ior and tunable materials properties (e.g., thermal, mechanical, optical, rheological) that 
are the basis of many industrial and technological applications. Thermally excited (Brow- 
nian) motion of nm-/^m-sized particles dispersed in a fluid medium drives the self-assembly 
of ordered phases. Important examples are nanoscale structures [4] and colloidal crys- 
tals, whose diverse crystalline symmetries and variable lattice constants can conveniently 
template photonic band-gap materials [5-7]. 

Interparticle interactions and correlations determine the distribution of microions 
(counterions and salt ions) around the colloidal macroions and thereby govern microion- 
mediated electrostatic interactions between macroions. As explained by the landmark 
theory of Derjaguin, Landau, Verwey, and Overbeek (DLVO) [8,9], repulsive interactions 
between like-charged macroions can stabilize a suspension against aggregation induced 
by van der Waals attractive interactions [10,11]. Equilibrium and dynamical properties 
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of charged colloids depend sensitively on electrostatic interactions, which can be widely 
tuned by adjusting system parameters: size, charge, and concentration of ion species; pH 
and dielectric constant of solvent. 

A powerful approach to modeling charged colloids is molecular simulation of either 
the primitive model, which includes all ions explicitly, or all-atom models, which include 
even individual solvent (e.g., water) molecules. At such microscopic resolution, simu- 
lations can potentially yield valuable insights into the thermodynamic, structural, and 
dynamical properties of real suspensions. Currently available processors and algorithms 
are limited, however, to relatively small systems and low ion size and charge ratios. Fur- 
thermore, simulations of microscopically detailed models do not necessarily elucidate 
the physical mechanisms underlying complex cooperative behavior. 

An alternative to brute-force modeling is simulation of a one-component model, de- 
rived by pre-averaging over the microion coordinates to obtain effective interactions be- 
tween macroions, screened by implicitly modeled microions. This coarse-grained strat- 
egy finesses the computational challenges that plague more explicit models, but relies on 
practical and accurate approximations for the effective interactions. A previous simula- 
tion study [12] validated the one-component model — implemented with linear-response 
and mean-field approximations — by direct comparison with pressure data from primi- 
tive model simulations at modest electrostatic couplings. The eventual breakdown of the 
model at higher couplings was attributed to failure of linear-screening approximations to 
account for nonlinear counterion association near macroions. 

A recently proposed charge renormalization theory of effective interactions in charged 
colloids [13] addresses limitations of linear-screening approximations by incorporating 
an effective (renormalized) macroion charge into the one-component model. Calcula- 
tions based on a variational approximation for the free energy indicated the potential 
of the theory to accurately predict thermodynamic properties of suspensions well into 
the nonlinear-screening parameter regime. The present paper describes complementary 
Monte Carlo simulations designed to test predictions for both thermodynamics (osmotic 
pressure) and structure (pair distribution function) of deionized suspensions of highly 
charged colloids. 

The remainder of the paper is organized as follows. Section |2| first reviews the prim- 
itive and effective one-component models of charged colloids. Section |3| outlines the 
charge renormalization theory, which predicts renormalized system parameters: effec- 
tive macroion charge, volume fraction, and screening constant. Section HI describes our 
Monte Carlo simulations, which take as input the renormalized effective interactions. 
Section |5| compares our results for the pressure and pair distribution function of deion- 
ized suspensions with corresponding data from simulations of the primitive model and 
with predictions of a variational free energy theory. Excellent agreement is obtained, with 
minimal computational effort, over ranges of macroion charge, volume fraction, and elec- 
trostatic coupling strength. Finally, Sec.|6|closes with a summary and conclusions. 
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2 Models 

The primitive model idealizes the colloidal macroions as charged hard spheres of monodis- 
perse radius a and bare valence Zq (charge — Zqc), the microions as point charges, and the 
solvent as a continuum of uniform dielectric constant e. We limit present consideration 
to monovalent microions, for which microion-microion correlations are relatively weak, 
and further assume index-matching of macroions and solvent, justifying neglect of image 
charges and polarization effects. At absolute temperature T, a characteristic length scale 
is the Bjerrum length Ag =e^/ (e/cgT), defined as the distance between a microion pair at 
which the electrostatic interaction energy rivals the average thermal energy. In a closed 
volume V, the suspension comprises Nm macroions and N± positive/ negative microions, 
of which Nc are counterions and Ng are salt ion pairs. Global electroneutrality relates the 
ion numbers via Nc = N+ — N_ = ZoNm- Alternatively, the suspension may be in Donnan 
equilibrium (e.g., across a semi-permeable membrane) with an electrolyte reservoir that 
fixes the salt chemical potential. 

While the primitive model has the virtue of explicitly representing all ions, simu- 
lations face severe computational challenges posed by long-range Coulomb interparti- 
cle interactions and large size and charge asymmetries between ion species. Sophisti- 
cated numerical algorithms, such as cluster moves [14, 15] and Ewald or Lekner sum- 
mation [16, 17], can significantly broaden the range of accessible length and time scales. 
Nevertheless, currently feasible studies are still limited to relatively small systems and 
modest asymmetries. 

An alternative approach to modeling charged colloids is derived from the Hamilto- 
nian H of the primitive model by formally tracing over the degrees of freedom of the 
microions (^) in the partition function of the multi-component mixture: 

Z=((exp(-HABT)),), = (exp(-HeffABT))^ , (2.1) 

leaving an explicit trace over only the macroion (ni) degrees of freedom. The resulting 
one-component model (OCM) is governed by an effective Hamiltonian 

Heff = E + x E + . (2-2) 

a function of the macroion coordinates r,, which comprises a one-body volume energy E, 
an effective macroion-macroion pair potential ^eff (^)' ^rid higher-order terms involving 
sums over effective many-body potentials. The volume energy, although independent 
of macroion coordinates, depends on the average densities of macroions and salt ions. 
The effective pair potential arises from screening of the bare Coulomb potential by the 
implicitly modeled microions. 

The computational benefit of reducing the number of components and the interac- 
tion range comes at a cost of increased complexity in the effective interactions. Linear- 
screening and mean-field approximations, which are justifiable in the case of weakly 
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Figure 1: Model of charged colloidal suspension: spherical macroions of radius a and point microions dispersed 
in a dielectric continuum. Strongly associated counterions in a spherical shell of thickness S renormalize the 
bare macroion valence Zq to an effective (lower) valence Z. 



correlated monovalent microions, yield analytical forms for the volume energy and ef- 
fective pair potential [18-20]. Moreover, many-body effective interactions are usually 
weak [21, 22]. The range of validity of the conventional OCM is nevertheless limited 
by the inevitable onset of nonlinear screening with strengthening macroion-counterion 
correlations. Previous studies [12,22] establish a parameter threshold ZqAb/^—'^ — char- 
acteristic of highly charged latex particles and ionic surfactant micelles — above which 
nonlinear effects tend to become significant. 

The OCM can be extended to more strongly correlated systems by incorporating 
physical concepts from charge renormalization theory [23]. Counterions that venture suf- 
ficiently close to a highly charged macroion may become thermally "bound" — closely 
associated with the macroion, if not condensed onto its surface [24,25]. As Fig. [l] depicts, 
counterions localized within a spherical shell surrounding the macroion may be consid- 
ered to renormalize the bare macroion valence. The composite "dressed" macroion, com- 
prising a bare macroion and its shell of bound counterions, carries a reduced effective 
valence Z < Zq. 

Counterions bind to a macroion at a distance at which the electrostatic energy of at- 
traction is comparable to the average thermal energy per counterion. Denoting by (p{r) 
the electrostatic potential at distance r from a macroion center, the width of the associa- 
tion shell 3 can be defined via 

e\cp{a+3)\ = CkBT , (2.3) 
where C is a dimensionless parameter of order unity. In general, the reference potential 
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should be chosen as the Donnan (average) potential (po of the suspension, C then being 
a function of salt concentration. Counterions within an association shell {a < r < a + 3) 
are presumed bound by the potential well of the respective macroion, while more distant 
counterions roam free and screen the dressed macroions. In other studies, a similar ther- 
mal criterion has been applied to the electrostatic potential [26-28] and the effective pair 
potential [29]. 



3 Theory 

The effective interaction theory proposed in ref. [13] combines a charge renormalization 
theory for the effective macroion charge with a mean-field, linear-response theory of the 
OCM [18-20]. Here we briefly outline the theory, whose predictions for effective interac- 
tions are subsequently used as input to our simulations (Sec.|4|). 

By invoking a mean-field random-phase approximation for the response functions 
(partial static structure factors) of a microion plasma, linear-response theory yields ana- 
lytical expressions for effective electrostatic interactions. This approach proves formally 
equivalent to density-functional [30-34], extended Debye-Hiickel [35], and distribution- 
function [36,37] formulations of effective interaction theory [4]. Linear-response theory 
is also equivalent to linearized Poisson-Boltzmann (PB) theory, based on a consistent lin- 
earization of the PB equation and expansion of the ideal-gas free energy functional to 
quadratic order in the microion density profiles [38]. 

Linearized theories predict the reduced electrostatic potential generated by a single 
bare macroion as 

tp{r) = ^cp{r) = -Zo\Bzr, ' r>a, (3.1) 

1 + Kfl r 

where j6 = 1 / A:^ T, k = i/47rAB («+ + «-) denotes the bare Debye screening constant, and 
n± = N± I [V'(l — ?/)] are the mean number densities of microions in the free volume, i.e., 
the total volume V reduced by the fraction r] occupied by the macroion hard cores. The 
constraint of global electroneutrality ensures that K depends implicitly on the average 
densities of both macroions and salt ions. 

Adapting Eq. (|3.1|) to the charge renormalization model (Fig. [T]), the potential around 
a dressed macroion, of effective valence Z and effective radius a^b, becomes 

where it = Y^47rAB(n_|_ + n_) now represents the effective (renormalized) screening con- 
stant, n± = N±/\V{l — fl)] and N± are the mean number densities and numbers of free 
microions, and fj = r]{\+5 / a)^ is the effective volume fraction of the dressed macroions. 
Combining Eqs. (|2.3P and (|3.2p yields the transcendental equation 

:C, (3.3) 
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which may be solved for the association shell thickness 5 for given values of Z,r], and C, 
noting that ii depends self-consistently on S. 

The distinction between free and bound microions implies a corresponding separa- 
tion of the total free energy 

^ — ^lree~^^o\md~^^mr (3-4) 

into three terms representing, respectively, contributions from free and bound microions 
and from macroion effective interactions. The free microions are, by construction, suf- 
ficiently weakly correlated with the macroions to be well described by linear-response 
theory. Their free energy (per macroion) is therefore accurately approximated by the 
linear-screening prediction for the volume energy (with renormalized parameters): 

w^=E/.K«.A')-il-f irl^-f WSt' 

where x± = N± / Nm are the free microion concentrations, A is the microion thermal de 
Broglie wavelength, and = Nm/V is the mean macroion number density. The three 
terms on the right side have physical interpretations as the ideal-gas (entropic) free en- 
ergy of the free microions, the self-energy of the dressed macroions, and the Donnan 
(average) potential energy of the microions. Although the bound counterions, being rela- 
tively strongly correlated with the macroions, do not yield to a linear-response treatment, 
their free energy (per macroion) can be reasonably approximated by 



^/bound = (Zq-Z) 



ln(^A3 



^^^^ 

+ (3.6) 



where the first term on the right side is an ideal-gas contribution, {Zo — Z)/Vs being the 
mean counterion density in the association shell of volume Vs = (47r/3) [{a-\-S)^ — a^], and 
the second term is the electrostatic energy, assuming tight localization of the counterions 
near r = a. 

The effective valence Z is now prescribed, for a given bare valence Zq, by the varia- 
tional ansatz [39-41] 

^ (/free + /bound) =0' (3-7) 



which ensures equality of the chemical potentials of free and bound counterions, under 
the constraint that Z and S are connected by Eq. ( I3.3D . The effective valence and shell 
thickness then determine the effective screening constant k. 

The macroion free energy F„, in Eq. I|3.4|) depends on correlations and effective inter- 
actions between dressed macroions. Incorporating renormalized system parameters into 
linear-response theory [19, 20], the effective pair potential is given by 

fiv,ii{r)=Z'\B[^:p^j r>2(fl + <5), (3.8) 
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whose screened-Coulomb form is identical to the long-range limit of the DLVO pair po- 
tential [9]. The renormalized effective interactions and system parameters can be input 
into computer simulations or liquid-state theories [42] of the OCM to compute thermo- 
djoiamic and structural properties of bulk suspensions of charged colloids. In the next 
section, we compare our simulation results with predictions of a variational approxima- 
tion [22, 30] for the macroion free energy (per macroion) based on first-order thermody- 
namic perturbation theory with a hard-sphere reference system: 

fm{nm,n±) =rnm\ fiis{nm,n±;d)+2nnrn [ drr^gyis{r,nm;d)vefi{r,nm,n±) \ . (3.9) 

(d) { Jd ) 

Here /hs and ^hs are, respectively, the excess free energy density and pair distribution 
function of the HS fluid, computed from the near-exact Carnahan-Starling and Verlet- 
Weis expressions [42]. Minimization with respect to the effective hard-sphere diameter d 
yields a least upper bound to the free energy. In practice, the renormalized system param- 
eters (Z, 3, k) must be held fixed in this minimization and in all partial thermodynamic 
derivatives. The corresponding prediction for the thermodynamic pressure is 

where the first three terms on the right side are ideal-gas contributions from macroions 
and free microions, the fourth term arises from density dependence of the self-energy of 
the dressed macroions, and the final term is generated by effective interactions between 
pairs of dressed macroions. While this simple variational theory yields the pressure and 
other thermodynamic properties, computer simulations or integral-equation theory are 
required to determine structural properties. 



4 Monte Carlo Simulations 



Working within the canonical (constant-NVT) ensemble, we consider a one-component 
fluid of macroions in a cubic box, subject to periodic boundary conditions, at fixed tem- 
perature, volume, macroion number, and mean salt concentration. The macroions in- 
teract with one another via effective electrostatic interactions that depend on the mean 
densities of both macroions and (implicitly modeled) salt. According to the standard 
Metropolis algorithm [16, 17], trial particle displacements are accepted with probability 

P = min{exp(-/3AU),1} , (4.1) 

where AU = U{n) — U{o) is the change in pair potential energy between the new (n) and 
old (o) states and 

Nn, 

U= £ V,ii{\Ti-Tj\) (4.2) 

i<i=l 
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is a sum of hard-sphere-repulsive- Yukawa (screened-Coulomb) pair potentials. Note that 
the acceptance probability for trial displacements does not involve the volume energy, 
since the mean density is fixed. Consequently, the volume energy does not affect the 
macroion structure, although it does contribute to the pressure (see below). To achieve 
high numerical precision, pair interactions were cut off at a distance Yc ~ 20/k, i.e., 20 
effective screening lengths. The cut-off radius determined the minimum box side length, 
L = IYc, necessary to avoid interactions of a particle with its own periodic images. For 
a given volume fraction, the requisite number of macroions was therefore prescribed as 
N„,~(3?//47r)(40/Kfl)3. 

We performed a series of simulations, using renormalized effective interaction pa- 
rameters (effective macroion valence, volume fraction, and screening constant), starting 
from face-centered cubic crystal configurations. Trial moves were executed by randomly 
displacing macroions with a step size adjusted to yield an acceptance rate of about 50%. 
Following an equilibration phase of 10^ cycles, statistics were accumulated for pressures 
and pair distribution functions over the next 10* cycles {Nm x 10* particle displacements). 
Test simulations for larger systems confirmed finite-size effects to be negligible. Rela- 
tively modest computing resources were required, with typical runs on a single Intel 
Xeon-HT processor lasting 30, 120, and 200 hours for = 2400, 4000, and 5800 particles, 
respectively. 

In the constant-Nyr ensemble, the total pressure may be expressed as 

/3p = n. + n,+n_-|^±g^ + /3pp,,, (4.3) 

where ppair is generated by effective interactions between pairs of macroions and corre- 
sponds to the last term on the right side of Eq. ( 13.101 ). The pair pressure is computed from 
the virial expression for a density-dependent pair potential [43]: 

where Vint is the internal virial, the volume derivative term accounts for the density de- 
pendence of the effective pair potential, angular brackets denote an ensemble average 
over configurations in the canonical ensemble, and ptail corrects for cutting off the long- 
range tail of the pair potential. The internal virial is given by 

Vint = ^ r, • f , = X] ( 1 + ) Veiiirij ) , (4.5) 
i=i i<j=i 

where f, = — J^j^iv'^ff{rij) is the effective force exerted on macroion / by all neighboring 
macroions within a sphere of radius r^. The second term on the right side of Eq. (|4.4|) is 
computed as the ensemble average of 



9 



Finally, neglecting pair correlations for r > Vc, the tail pressure is approximated by 

Ptaii = - -J^mj^ drr^v^^^{r) = —n„, ( ^ 1 r^z^eff (f.) ■ (4-7) 

The macroion structure of the suspension is characterized by the macroion-macroion 
pair distribution function g{r), defined such that 4nr^g{r)dr equals the average number 
of macroions in a spherical shell of radius r and thickness dr centered on a macroion. 
With each particle regarded in turn as the central particle in a given configuration, neigh- 
boring particles were assigned to concentric spherical shells (bins), of thickness Ar = 0.1a, 
according to their radial distance r from the central particle. Following equilibration, 
g{r) was computed in the range 2{a+S) < r < L/2 by accumulating statistics of particle 
numbers in radial bins and averaging over all configurations. The raw distributions were 
finally smoothed by averaging each bin together with its neighbors in a moving average 
algorithm. 



5 Results and Discussion 

The renormalized effective interactions (Sec. |3]) previously provided the basis for vari- 
ational theory calculations of the pressure of deionized suspensions [13]. The same 
effective interactions are here input into Monte Carlo simulations (Sec. H]) of the one- 
component model to compute both the pressure and macroion-macroion pair distribu- 
tion function. Comparing our results with corresponding data from primitive model 
simulations allows testing the accuracy of the effective interactions and the variational 
free energy approximation. 

All results presented below are for the case of monovalent counterions, zero salt con- 
centration, and an aqueous solvent at room temperature (Ag =0.714 nm). Throughout, the 
dimensionless thermal parameter in Eq. (|3.3|) is fixed at C = 3, a value shown in ref . [13] 
to give satisfactory agreement with pressure data from primitive model simulations. For 
salt-free suspensions, this choice of C corresponds toe|^(fl+^) — ^d| = 2A:b T. In passing, 
we note that the parameter C is analogous to the adjustable cell radius b in the PB cell 
theory of Zoetekouw and van Roij [44]. 

Figure |2] illustrates the distinction between the bare and effective macroion valences. 
For a sufficiently small bare valence, Eq. (|3.3|) has no real-valued solution for any nonzero 
association shell thickness. In this case, all counterions are free (5 = 0, Vs = 0) and the free 
energy is minimized by Z=Zo. Beyond a threshold bare valence, however, the association 
shell appears and rapidly thickens with increasing Zq. Correspondingly, the free energy 
minimum shifts to Z < Zq . With increasing Zq, the effective valence grows logarithmically, 
in contrast to the saturation observed for polyelectrolytes [24, 25] and predicted for col- 
loidal suspensions by PB cell theories [23,45,46] and Debye-Hiickel theories [39-41,47,48]. 

We first test the capacity of the charge renormalization theory [13] to predict the 
macroion structure. Figures |3] and H] compare results for the pair distribution function 
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Figure 2: Effective valence Z vs. bare valence Zg for a deionized suspension (cs~0) of macroions of radius a = 2 
nm and volume fraction ;^ = 0.01. Inset: counterion association shell emerges and thickens beyond threshold Zq. 



g{r) from our simulations of the OCM, using charge-renormalized effective interactions 
as input, and corresponding data from extensive Monte Carlo simulations of the primi- 
tive model (PM) [49] for salt-free suspensions over ranges of bare valence, volume frac- 
tion, and electrostatic coupling strength F = Ag/a. Well beyond the threshold for charge 
renormalization {ZqY ~ 7), the OCM and PM results agree closely. Significant deviations 
emerge only for ZqT > 15, the OCM consistently overpredicting the macroion structure. 
For reference. Fig. |3] also includes OCM results obtained with unrenormalized effective 
interactions, illustrating the impact of charge renormalization on strongly coupled sus- 
pensions. 

Finally, we test predictions of the charge renormalization theory [13] for thermo- 
dynamics. Figures |5] and |6| directly compare the pressures computed from our OCM 
simulations [via Eqs. (|4.3 |) - (|4.7|) 1 with corresponding data from primitive model simu- 
lations [49] of salt-free suspensions over ranges of bare valence, volume fraction, and 
electrostatic coupling strength. The OCM and primitive model results are seen to be in 
excellent agreement up to the highest couplings for which primitive model data are avail- 
able (ZoF~30), well beyond the charge renormalization threshold. Figures |5| and |6| also 
compare the pressures resulting from our OCM simulations with predictions of the vari- 
ational theory [Eq. I|3.10|l l. The near exact agreement between simulation and theory for 
the OCM confirms the accuracy of the variational approximation for the macroion free 
energy [Eq. (|3.9|) 1 and demonstrates the equivalence of the thermodynamic and virial 
expressions for the pressure, i.e., Eqs. (|3.10|) and (|4.3|) . respectively. 
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Figure 3: Macroion-macroion pair distribution function g{r) vs. radial distance r (units of macroion radius a) 
of salt-free suspensions computed from Monte Carlo simulations of the effective one-component model, with 
(0CM2, dashed curves) and without (OCMl, dotted curves) charge renormalization, and of the primitive model 
(PM) [49] with explicit counterions (solid curves) for volume fraction rj = O.Q\, electrostatic coupling parameters 
I=XB/a= 0.0222, 0.0445, 0.0889, 0.1779, 0.3558, 0.7115 (bottom to top), and bare macroion valence Zo=40 
(a), 80 (b), and 160 (c). For clarity, curves are vertically offset in steps of 0.5. 
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Figure 4: Macroion-macroion pair distribution function g{r) vs. radial distance r (units of macroion radius 
a) of salt-free suspensions computed from Monte Carlo simulations of the charge-renormalized effective one- 
component model (0CM2, dashed curves) and of the primitive model [49] (PM, solid curves) for various bare 
macroion valences Zq, electrostatic coupling parameters T = \g/a, and volume fractions rj= 0.00125, 0.0025, 
0.005, 0.01, 0.02, 0.04, and 0.08 (right to left). 
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Figure 5: Total reduced pressure /3p/ntot vs. macroion volume fraction rj, where Mtot = (Zo + l)wm (total 
ion density), of salt-free suspensions with bare macroion valence Zo=40 and electrostatic coupling constants 
(from top to bottom) r= 0.0222, 0.0445, 0.0889, 0.1779, 0.3558, 0.7115. Open symbols: Monte Carlo 
simulations of the primitive model [49]. Filled symbols: Monte Carlo simulations of the effective one-component 
model. (Symbol sizes exceed error bars.) Curves: variational theory with (solid) and without (dashed) charge 
renormalization. The double-ended arrow points to corresponding curves for r = 0.3558. The dashed curve for 
r = 0.7115 is off-scale, the pressure being negative. 
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Figure 6: Total reduced pressure /Jp/wtot vs. electrostatic coupling parameter T of salt-free suspensions with 
fixed volume fraction r] = 0.01 and bare macroion valence (top to bottom) Zq = 10, 20, 40, 80. Open symbols: 
Monte Carlo simulations of the primitive model [49]. Filled symbols: Monte Carlo simulations of the effective 
one-component model. (Symbol sizes exceed error bars.) Curves: charge-renormalized variational theory. Inset: 
Ratio of effective to bare macroion valence Z/Zq vs. T. 
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6 Conclusions 

Summarizing, we have performed Monte Carlo simulations of the one-component model 
of charged colloids, using as input effective electrostatic interactions predicted by a re- 
cently proposed charge renormalization theory [13]. Structural and thermodynamic prop- 
erties of salt-free suspensions are computed and directly compared with simulations of 
the primitive model. For bare valences Zq and electrostatic coupling strengths T consid- 
erably exceeding the renormalization threshold, the pair distribution function and pres- 
sure are in close agreement with corresponding results from simulations of the primitive 
model. Significant discrepancies between the OCM and PM simulations are observed 
in g{r) for ZqT > 15, and in the pressure only for ZqT > 30. The level of agreement ap- 
pears to be comparable to the charge-renormalization scheme of ref. [14], which is based 
on a PB cell model [23]. The computationally practical approach described here may 
provide a useful modeling alternative, for some applications, to molecular-scale simu- 
lations of charged colloids. Further comparisons with primitive model simulation data, 
particularly at nonzero salt concentrations, would help to chart the validity range of the 
charge-renormalized OCM. 

The threshold for charge renormalization closely coincides with the onset of a spin- 
odal phase instability predicted by linearized effective-interaction theories [22, 30, 32, 35] 
for deionized (but not salt-free) suspensions of highly charged colloids. Although obser- 
vations consistent with bulk phase separation have been reported [50-57], the experimen- 
tal picture of deionized suspensions remains cloudy. Furthermore, the relevant param- 
eter regime is not yet accessible to primitive model simulations. Recent studies based 
on the PB cell model have concluded that the predicted instability may be an artifact 
of linearization approximations [58-60]. Preliminary calculations based on the present 
approach indicate, however, that the phase instability may survive charge renormaliza- 
tion [61]. Future work will continue to explore the remarkable phase behavior of deion- 
ized suspensions. 
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